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We study the problem of efficient integration of variational equations in multi-dimensional 
Hamiltonian systems. For this purpose, we consider a Runge-Kutta-type integrator, a Taylor 
series expansion method and the so-called 'Tangent Map' (TM) technique based on symplectic 
integration schemes, and apply them to the Fermi-Pasta-Ulam /3 (FPU-/3) lattice of N nonlin- 
early coupled oscillators, with N ranging from 4 to 20. The fast and accurate reproduction of 
well-known behaviors of the Generalized Alignment Index (GALI) chaos detection technique is 
used as an indicator for the efficiency of the tested integration schemes. Implementing the TM 
technique-which shows the best performance among the tested algorithms-and exploiting the 
advantages of the GALI method, we successfully trace the location of low-dimensional tori. 

Keywords: Hamiltonian systems, numerical integration, variational equations. Tangent Map 
method, GALI method 



1. Introduction 

From interactions of stars in galaxies to particle beams in high energy accelerators, Hamiltonian mechanics 
is found at the very heart of modehng and understanding dynamical processes. The necessity to classify 
evermore complex systems in terms of stabihty and pred ictability has lead to a wealth of methods discrim- 
inating chaotic from regular behavior (see for example Skokod . 2010l . Sect. 7]). Most of these techniques 
rely on the study of the time evolution of deviation vectors of a given orbit to discriminate between the 
two reg imes. The time evolution of these vectors is governed by the so-called variational equations. 

In (Skokos fc Gerlachl . bold ] the 'Tangent Map' (TM) technique, an efficient and easy to implement 
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method for the integration of the vari ational equations o f Harniltonian systems based o n the use of sym- 
plectic integrators was introduced. In [Skokos fc Gerlachl . l20ld : iGerlach fc Skokod . l201(]| | the TM method 
was apphed mainly to low-dimensional Hamiltonian systems of 2 and 3 degrees of freedom, and proved to 
be very efficient and superior to other commonly used numerical schemes, both with respect to its accuracy 
and its speed. 

The scope of the present work is to extend these results by investigating whether the efficiency of the 
TM method persists also when multi-dimensional Hamiltonian systems are considered. The study of such 
systems presents a challenging numerical task, which makes the use of fast and accurate numerical tools 
imperative. In the present paper we use as a toy model the famous Fermi-Pasta-Ulam /3 (FPU-/3) lattice, 
which is presented in Sect.O In Sect.[3]the different numerical methods we use, i.e. the Generalized Align- 
ment Index (GALI) chaos indicator, the TM method and the SABA family of symplectic integrators, the 
Taylor series integrator TIDES, and the general-purpose high-accuracy Runge-Kutta integrator DOP853, 
are presented and their properties are briefly discussed. Then, in Sect. 14.11 the numerical results of the 
application of these numerical schemes for the integration of variational equations of the FPU system are 
presented, while in Sect. I4.2l the GALI method is used for locating motion on low-dimensional tori. Finally, 
in Sect. Owe summarize our results. 



2. The FPU lattice 



As a mode l of a multi-dimensional H amiltonian system we consider the FPU-/? lattice [Fermi et all Il955l : 
Ford . 1992 : Chaos Focus Issu3 . 2005], which describes a chain of particles with nearest neighbor inter- 
action. Regarding numerical integration algorithms, the FPU lattice is a very challenging problem, since it 
exhibits oscillations on largely different time scales. The Hamiltonian of this degrees of freedom (AD) 



system as a function of the momenta p 



Hn = HN{p,q) 



{pi,...,PN 

N 



and the coordinates q = (gi, . . . , qn) is given by 
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In our study we impose fixed boundary conditions, 



1. e. 



+ 



qN+i = 0, set /3 = 1.5, and consider 



models whose number o f particles vary from N 
Skokos fc Gerlachl . l20inl | the particular case A^ 



4 up to A^ = 20. We note that in [Skokos et al. 
> was studied in detail. 



2008 



The Hamiltonian ([T]) can be split into two parts A and -B, which respectively depend only on the 
momenta and the coordinates, i. e. H]\j = A{p) + B{q). The Hamilton's equations of motion are 



Pi 



and 



Pi 



opi oqi 

with 1 < i < N , and 5j denoting the Kronecker delta, which is equal to 1 if i = j and to otherwise. The 
variational equations governing the time evolution of a deviation vector w = {5qi, . . . , 6qiy, 6pi, ■ ■ ■ ,5pn), 
that evolves in the tangent space of the Hamiltonian's phase space are given by 

Spi = Sqi and 6qi = - 'S'^ 5qj . (3) 



dqidqj 



3. Numerical methods 

In order to compute the stability of a particular solution of Hamiltonian ([T]) , or in other words of an orbit in 
the 2A^-dimensional phase space of the system, the equations of motion Q have to be integrated together 
with the variational equations The time evolution of the latter contains information on the stability 
of the orbit, which can be quantified by using some chaos indicator, for example the GALIs. Different 
numerical approaches can be used to solve the system of ordinary differential equations given by Eqs. ^ 
and ([3]). In this section, after briefly recalling the definition of the GALI and its behavior for regular 
and chaotic motion, we present an overview of the methods used in the curre nt study. A more detailed 
description of further possibilities based on symplectic methods can be found in [Skokos fc Gerlachl . l20inl |. 
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3.1. The Generalized Alignment Index (GALI) 

The GALI was originally introduced in [Skokos et alV |2007| | as an efficient c haos de t ection method, gen- 



eraliz i ng a similar indicator called the Smaller Alignment Index (SALI) (Skokod . 2001 : Skokos et al. 



2003I . 2004 1 ■ The method has been applied successfully to different dynamical systems for the discrimi 

r the detection of regular motion on low dimen- 
20081 : iBountis e,t oil . l2009l : iManos fc RiTffol . I2OI0I : 



sional tori Christodoulidi & 


Bountis. 2006: Skokos et al. 


Manos k, Athanassoulal. 2011 


: Manos et al. 


. 20111. 



For A'^D Hamiltonians the Generalized Alignment Index of order k (GALI^), 2 < k < 2N, is determined 
through the evolution of k initially linearly independent deviation vect o rs Wk (0), which are continually 
normalized, keeping their directions intact. According to Skokos et al\ [2007l | GALI^ is defined as the 
volume of the fc-parallelepiped having the k unitary deviation vectors Wi{t) = Wi{t)/\\wi{t)\\, i = 1,2, . . . , k, 
as edges. GALI^ is therefore determined through the wedge product of these vectors 



GALIfc(t) = A W2it) A • • • A Wkit)\\, 



(4) 



with II • II denoting the usual norm. The behavior of GALI^, for regular motion on a n s-dimensional torus 
is given by [Skokos et all . l2007l : IChristodoulidi k Bountii . I2OO6I : ISkokos et all . l2008l | 



GALIfc(t) oc < 



' constant if 2 < k < s 

iis <k<2N -s 
if2N-s<k<2N 



(5) 



while for chaotic orbits GALI^, tends to zero exponentially following the law [Skokos et all |2007| | 

GALIfc(t) oc e-[('^i-'^2)+('^i-'^=^)+-+('^i-'^'=)l*, (6) 

where ai, . . . ,ak are the first k largest Lyapunov characteristic exponents of the orbit. 

We chose to apply the GALI method in order to check the accuracy of the numerical integration of 
variational equations because this index depends on the evolution of an ensemble of deviation vectors. Thus, 
even small errors in the integration of these vectors are expected to significantly influence the evolution 
of GALIs. In particular, we focus our attention on regular orbits lying on tori of various dimensions s, for 
which the GALIs exhibit different evolutions depending on s and the order k of the index. The possible 
deviations of numerically evaluated GALIs from their expected behaviors ([5]) will identify the inaccurate 
integration of the deviation vectors far better than it would be the case using chaotic orbits, for which all 
GALIs tend to zero exponentially ([6]). Throughout the paper we denote regular orbits on an s-dimensional 
torus of the A^D system ([T]) as 7^^, where s can range from 2 to N. 



3.2. The Tangent Map method using symplectic algorithms 

Symplectic methods are often the preferred choice when integrating dynamical problems, which can be 
described by Hamiltonian functions. A thorough discussion of such methods can be found in Hairer et al\ 
(2002l |. Let us just mention some properties of symplectic integrators which are of interest for our study. 
Symplectic methods cannot be used with a trivial automated step size control. As a consequence, they 
are usually implemented with a fixed integration step r. Due to their special structure they preserve the 
symplectic nature of Hamilton's equations intrinsically, which in turn leads to results that are more robust 
for long integration times. A side-effect of structure preservation is that the error in energy remains bounded 
irrespe ctive of the total integrat ion time. 

In Skokos &: Gerlachl . 12010 ] it was shown that it is possible to integrate the Hamilton's equations of 
motion and the corresponding variational equations using the TM technique based on symplectic splitting 
methods. Let u s outline the basic idea be hind the TM method, which is founded on a general result stated 
for example in Laskar &: Robutel . 2001 1: Symplectic integrators can be applied to systems of first order 
differential equations X = LX, that can be written in the form X = [La + Lb)X, where L,La,Lb are 
differential operators defined as L-^f = {x, f} and for which the two systems X = LaX and X = LbX 
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are integrable. Here {f,g} are Poisson brackets of functions f{q,p), g{q,p) defined as: 

{f.g} 



N 

E 

1=1 



df dg 
dpi dqi 



df dg 
dqi dpi 



(7) 



The set of Eqs. ([2]) is one example of such a system, since the Hamiltonian ([T]) can be divided into two 
integrable parts A and B with H = A(p) + B{q) as already noted. A symplectic splitting method separates 
the equations of motion ([2]) into several parts, applying either the operator La or Lb- These are the 
equations of motion of the Hamiltonians A and B, which can be solved analytically, giving explicit mappings 
over the time step qt, where the constants Cj are chosen to optimize the accurac y of the integrator. These 
mappings can then be combined to approximate the solution for a time step r. In Skokos &: Gerlachl . l201Cll | 
it was shown that the derivatives of these mappings with respect to the coordinates and momenta of the 
system (the so-called tangent maps) can be used to calculate the time evolution of deviation vectors or, in 
other words, solve the var i ationa l equations ([3|). 

In (Laskar &: Robutel . 200 1| a family of symplectic splitting methods called SABA„ and SBAB„ was 
introduced, with n being the number of applications of operators La and Lb- These integrators were 
designed to have only positive int ermedia t e step s. Since it is not possible to construct symplectic integrators 
of ordei0 > 2 with this property Suzukil . 1991], small negative corrector steps C can be added before and 
after the main body of the integrator to further increase the accuracy. In Sect. 14.11 we test 3 different 
integrators, namely SABA2, SABA2C and SBAB2C, present the results of this comparison and discuss the 
characteristics of these algorithms when applied to the FPU system. 

The fact, that standard adaptive stepping policy is not possible with symplectic integration schemes 
necessitates an initial assessment of s tability for the algo rithms used, in order to derive an upper bound 
on the choice of step-size. Following iHairer et al\ 2002l | we note that SBABi is equivalent to the well 
known St0rmer-Verlet or leap-frog method, and SABAi to its adjoint. In order to perform a linear stability 
analysis of the FPU-/3 lattice pro blem for SABA2, we introduce normal mode momenta Pi and coordinates 
Qi as it was done for example in [Skokos et al . 2008]. The unperturbed (/3 = 0) Hamiltonian ([1]) can then 
be written as a sum of the so-called harmonic energies Ei, i.e. 

TV N 



2sin 



ITT 



2{N + 1] 



(8) 



where are the corresponding harmonic frequencies. With M denoting the Jacobian of the numerical 
mapping from initial coordinates and momenta to those at the next step for SABA2 we have 



M 



Qiir) 



M 



Qi{0) 



1 



8V3 



The characteristic polynomial of matrix M amounts to A + A I —2 + 



2v^ 



24 



the eigenvalues to be of modulus one, given the harmonic frequencies dS]) satisfy ]ci;i] < 2, 



the maximum 



admissible step-size for SABA2 will be Tn 



1.6. Thus, in our study we always use r < Tn 



3.3. Taylor series methods -TIDES 

The basic idea of the so-called Taylor series methods is to approximate the solution at time + r of a 
given d-dimensional initial value problem 

dX{t) 



dt 



f{X{t)) XGM^^G 



(9) 



^In this work we call a symplectic integrator to be of order n, if it introduces an error of 0{t") in the approximation of the 
real solution, with t being the integration time step. 
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from the n^^ degree Taylor series of X{t) at t = tf. 



X{ti+T) ~ X{ti)+T 



+ 



d^XiU 



+ . . . + 



(10) 



dt '2! dt2 n\ dt" 

For details see e.g. Hairer et all 19931 . Sect. 1.8] and references therein. In this work we call an integrator 
being of order n, when the first neglected term in this Taylor series expansion is of 0{t"'^^). The compu- 
tation of the derivatives can be very cumbersome, depending on th e struc ture of /, and is done efficiently 
using a utomatic differentiation techniques (see for example [B arriol . 2005]). 

In Gerlach &: Skokoi . 2010( ] two different publicly available implementations of the Taylor method were 
used and compared regarding thei r reliab i lity and efficiency i n the case of a 2D Hamiltonian system. One of 
these integrators, called TIDE^ Barrio . 2005 : Abad et al. , 2010l ]. which showed better performance, will 
also be used in this work as a representative of methods based on Taylor series expansions. TIDES comes 
as a Mathematica notebook. After inserting the differential equations one desires to integrate, the notebook 
generates automatically all the necessary subroutines to compute the given problem. The FORTRAN code 
produced by TIDES was included into the existing testbed of our work without further modification. 

In order to obtain optimal results, the TIDES algorithm is free to choose its order and step size during 
the whole integration interval. We used one parameter only, the so-called one step accuracy 6, to control 
the numerical performance of the algorithm. To be more precise, an integration step from time ti to tj + r is 
accepted, if the internal accuracy checks estimate that the local truncation error of the solution X(ti + r) 
is less than 6. If this error is too large, the integrator automatically tries to increase the internal order 
and/or adjust the step size r. After each successful step the deviation vectors are re-normalized. 

Let us remark that another elegant way to express Eg. (1 10 II can b e achieved using the so called Lie 



series formalism. Lie series ha ve been rediscovered by [Grobneri. 



of dynamical astronomy (e.g. Hanslmeier &: Dvorakl . 19841 : Delva . 1 1983 ]) to numerically solve ordinary 



967| and used extensively in the field 



differential equations. Sharing the same ansatz with symplectic maps (Sect. 13. 2p . Lie series can be used to 
iterate first order ordinary differential equations of the form X = LX as follows: 



X 



t+T 



e^Xt 



j=0 



Xi + 0(r"+i) 



fill 



Note that contrary to Sect. 13.21 no assumptions on the properties of the differential operator are made. 
Thus, by evaluating the consecutive derivatives LX, L'^X, Li^X and building the corresponding exponen- 
tial series up to order n one is able to follow the trajectory of X through phase space. The truncated Lie 
series' numerical map Xt Xt+r is not area-preserving, in general, since the method is non-symplectic. 
Therefore, the implementation of an adaptive step-size into Lie series algorithms is possible without reser- 
vations. This can be achieved v ia estimates on the size of the remainder of the exponential series (see for 
example Eggl fc Dvorakl . I20 id ]). 

Algorithms using series expansions become most efficient when combined with recursion relations, 
where higher order derivatives can be calculated from previous ones. In the course of this work, also 
extensive tests have been undertaken to adapt the Lie series method to the FPU-/3 lattice. Since for 
this problem such recursions are not available, we used algebraic manipulation software to compute the 
successive applications of the operator L to Xt, and implemented the generated code into a FORTRAN 
program. Since this approach proved to be reliable but computationally very expensive, we do not include 
it in the discussion of our results in Sect. HI 



3.4. General purpose integrators-DOP853 

In general, the computation of higher derivatives of functions X{t) soon becomes very complicated. There- 
fore, the methods described in Sect. I3.3l became popular only recently, after automatic differentiation could 
be performed efficiently by computers. Before that, other methods were developed to approximate Eq. (jlOp . 



^Freely available at http://gme.uiiizar.es/software/tides. 
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One of these are the so called Runge-Kutta methods (see for example jHairer et al . 1993 . Sect. II. 1.1]). 
An s-stage Runge-Kutta method is given as 

s 

X{ti + T) = X{ti) + T'^biki with (12) 

i=l 

s s 

h = f{ti + CiT,X{ti) + r'^aijkj) and Cj = ^ay. 

j=i i=i 



The real numbers 6j, Ojj with i, j = 1, . . . , s are chosen to approximate Eq. (flOj) to the desired order. If one 
requires further aij = for i < j the integration method will be explicit. Runge-Kutta integrators exist as 
symplectic as well as non-symplectic variants. 

In this work a 12-stage explicit Runge-Kutta integration method c alled DQP853 is us e This non- 
symplectic scheme is based on the method of Dormand and Price (see Hairer et al. . 19931 . Sect. II. 5] for 



further details). With this integrator we solve the set of differential equations composed of Eqs. ([2]) and 
(l3|) simultaneously. Here we use again the parameter 5 to control the integrator's overall behavior. For 
the DOP853 integrator the estimation of the local error and the step size control is based on embedded 
formulas of orders 5 and 3. 



4. Results 

Before investigating how variational equations can be integrated efficiently, one should first clarify the 
meaning of the term 'efficiency'. On the one hand, an efficient integration is one that is performed as fast 
as possible, and on the other hand, this computation should be also done as accurately as possible. Since 
accuracy always comes at the cost of intense computational efforts -meaning large CPU times- it tends 
to contradict the first mentioned aspect of efficiency. In addition, since we are especially interested in the 
computation of chaos indicators, an accurate computation should imply the correct distinction between 
regular and chaotic orbits of the studied dynamical system. Thus, in this work we consider the integration 
of variational equations to be efficient, when the computed GALIs are obtained with the least-possible 
CPU time requirements, achieving at the same time the correct characterization of orbits as regular or 
chaotic. 

In the next section we present a thorough discussion of how variational equations can be integrated 
numerically using the methods described in Sect. [31 We explain possibilities of estimating the accuracy of 
such computations and discuss the obtained results with respect to our definition of efficiency. In Sect. 14.21 
we apply this knowledge for a more global investigation of the properties of the FPU-/3 lattice. As a final 
remark we note that all presented computations were performed on an Intel Xeon X3470 with 2.93 GHz, 
using extended precision (80 bit). 



4.1. Efficient integration of variational equations 

In most applications it is common to integrate the variational equations together with the equations of 
motion they are based on. Thus, an important prerequisite to obtain correct stability results is the correct 
computation of the dynamics of the system itself. For this reason we first discuss the properties of the 
integrators described in Sect. [3] when applied to some specific orbits of the FPU-/? lattice. In general, 
an accuracy estimate of a numerically obtained orbit of a dynamical system can be given via monitoring 
how well the system's ffi'st integrals (e.g. the total energy, total linear momentum etc.) are conserved. 
For conservative Hamiltonian systems, like this can be done easily by checking the conservation of 
the energy H itself. The absolute value of the relative errors of the total energy \/S.H/H\, for different 
integrators and step sizes r, for a T2 orbit over t = 10^ time units are given in Tabled! For non-symplectic 
methods also the one-step accuracy 5 is mentioned. 



■^DOPSSS is freely available from http://www.unige.ch/~hairer/software.html. 
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Table 1. Information on the performance of the different numerical methods used 
for the computation of all the GALIs of the T2 orbit with initial condition Qi = 0.1, 
Pi = 0, 1 < i < 4, of system ((T| with A*' = 4, over t — 10^ time units. The given r for 
non-symplectic methods is computed as a mean step size r = t/ua, where Ua is the 
number of accepted steps, while 5 is the one-step accuracy. The absolute value of the 
relative energy error \AH/H\ at the end of each integration, as well as the required 
CPU time for each method are also reported. The last column provides information 
on whether the computed GALIs identified the nature of the regular orbit correctly 
(Y) or not (N) within the time interval t = 10^ (see also Fig. [T]). 



Integrator S r Order \AH/H\ CPU time Correctness 



TM-SABA2 






1.00 


2 


6 


X 


10" 


-2 





min 


03 sec 


Y 


TM-SABA2 






0.50 


2 


2 


X 


10" 


-3 





min 


06 sec 


Y 


TM-SABA2C 






0.50 


4 


4 


X 


10" 


-5 





min 


09 sec 


Y 


TM-SABA2C 






0.10 


4 


5 


X 


10" 


-7 





min 


46 sec 


Y 


TIDES 




5 


0.66 


10 


4 


X 


10" 


-2 





min 


45 sec 


N 


TIDES 


10" 


8 


0.60 


14 


1 


X 


10" 


-5 


1 


min 


14 sec 


N 


TIDES 


10" 


10 


0.54 


16 


2 


X 


10" 


-7 


1 


min 


37 sec 


Y 


DOP853 


10" 


5 


1.32 


8 


5 


X 


10" 


-1 





min 


11 sec 


N 


DOP853 


10" 


10 


0.25 


8 


8 


X 


10" 


-6 





min 


54 sec 


Y 


DOP853 


10" 


11 


0.19 


8 


6 


X 


10" 


-7 


1 


min 


11 sec 


Y 



Firstly, we focus on the application of the TM method with the SABA and SBAB integrators. The 
corresponding results are given in the first 5 lines in Table [TJ Comparing the energy conservation between 
SABA2 and SABA2C one finds that the use of the corrector steps significantly improves \AH/H\ for the 
same step size of r = 0.5. This result can of course be explained by the fact that SABA2 is an integration 
scheme of order 2, while SABA2C is of order 4. Reducing the step size for SABA2C to r = 0.1 further 
improves the energy conservation as expected. We note that this reduction leads to a linear growth by a 
factor 5 of the required CPU time, as expected. Since SABA2C shows the best performance, we use only 
this integrator for further investigations. 

Comparing the results of SABA2C with the ones obtained by the non-symplectic methods TIDES and 
DOP853 one finds that both non-symplectic methods need more CPU time in order to reach the same final 
value of \AH/H\. If one computes a mean step size for these algorithms, defined as the total integration 
time t{= 10^) divided by the number of accepted steps Ha, one finds that both integrators achieve this 
final energy error by a larger mean step size, compared to SABA2C, which is due to the higher orders of 
these integrators. While the highest order for DOP853 is fixed to 8, TIDES uses automatic order selection, 
which explains the larger mean time step for the same value of one-step accuracy 5. 

While monitoring energy conservation serves as a control parameter over the state vector of the system 
itself, it lacks information on how accurately the corresponding variational equations are solved. If the 
stability of certain initial conditions is known, one can use the theoretically predicted behaviors ([5]) of 
the GALI chaos indicator to estimate the reliability of the numerical computation. Therefore, in the last 
column of Table [H we provide information on whether the integration was able to identify the 7'2 orbit 
correctly as being a regular orbit lying on a 2-dimensional torus. 

The time evolution of GALI^, k = 2, ... ,8, for some of the runs of Table [1] is shown in Fig. [H From 
Eq. ([5]) it is known that GALI2 should be constant for a T2 orbit, while GALI3 and GALI4 should decrease 
proportionally to and respectively. A correct characterization is possible by using the TM method 
with SABA2 and a step-size as large as r = 1.0, although the corresponding energy error \AH/H\ ~ 10~^, 
is rather high (see first line in Table [1]). For the non-symplectic methods it is found that 6 = leads 
to a similar error in energy conservation, but is not sufficient for the correct computation of the GALIs 
(see the first and third columns of Fig. [T]) . Decreasing the one-step accuracy 5 improves the accuracy, but 
the integrations become less and less efficient compared to the TM integration as the required CPU time 
grows. We note that a correct dynamical characterization of the orbit is possible both for the DOP853 and 
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DOP853, 5=10"^ DOP853, 5=10"''° TIDES, 5=10"^ TIDES, 5=10"'' 




12345 12345 12345 12345 

lOQiot lOQiot lOQiot lOQiot 



Fig. 1. Time evolution of GALIs for the regular orbit with initial condition q.i — 0.1, Pi = 0, 1 < i < 4 of system ([T]) 
with A*' = 4, as computed using non-symplectic schemes. The results are given as colored curves. The TM method results with 
SABA2C and T = 0.5 are presented by grey-scale curves in the background serving as a reference. These curves are not always 
clearly visible as they are overlapped by the colored ones. 



TIDES for 6 < 10-^°. 

Using energy conservation as an indicator for the quality of an integration one could argue that the 
difference in CPU time between SABA2C with r = 0.1 and DOP853 with 5 = 10"^° is not very significant. 
While this is true for the T2 orbit, the difference becomes more pronounced, when the number of particles 
of the FPU-/3 lattice is increased. This is evident in Fig. [2] where we plot, as function of N, the ratio of 
the required CPU times between the TM method with SABA2C and the TIDES (blue curves) and DOP853 
(red curves) integrators for a regular orbit with initial condition qi = 0.1, Pi = 0, I < i < N for = 4, 

8, 12 and 20. If after t = 10^ time units an error in energy conservation of \AH/H\ ~ 10~^ is sufficient, 
one could use either SABA2C with r = 0.5, TIDES with 5 = 10"^ (although for A = 4 the corresponding 
GALIs do not exhibit the theoretically expected behavior for the whole time interval), or DOP853 with 
6 = 10"^'' (see Table [T]). The ratio of the CPU time of these runs is given in Fig. [2]^a). In Fig. EJ^b) we show 
a similar comparison when SABA2C with r = 0.1, TIDES with S = 10"^°, and DOP853 with 6 = 10"^^ 
are used, which yields \AH/H\ < 10~^ at the end of the integration. In this case, the GALIs computed by 
these methods show the time evolution predicted in Eq. ([5|). 

Figure [2] shows that the efficiency of the TM method improves with increasing N when compared to 
the non-symplectic methods used in this study. While the ratio is 2 between the CPU times of SABA2C and 
TIDES for A = 4 in Fig. [2]^b) it increases up to 7 when using A = 20 particles in the FPU-/3 lattice. The 
ratios are generally smaller when comparing the TM method with DOP853, but also here a growing trend 
can be observed. The results become even more pronounced when a larger error in energy conservation is 
acceptable, as it is shown in Fig. [2]J^a). 

When analyzing dynamical systems one very often has to follow a trajectory for long time intervals 
to determine the behavior of a specific initial condition correctly. Especially for weakly chaotic orbits 
differences to regular orbits are shown only in the late evolution of chaos indicators. In such investigations 
the TM method also proves to be superior compared to other techniques. As an example we show in Fig. [3] 
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Fig. 2. The ratio of required CPU time between integrations using the TM method with SABA2C and the TIDES (blue 
curves) and DOP853 (red curves) integrators for a TJ^2 regular orbit of Hamiltonian (ITJ with initial condition g,j — 0.1, Pi = 0, 
1 < i < TV for A'^ = 4, 8, 12 and 20. Comparisons are performed between (a) TM method with step size r — 0.5, TIDES 
with S = 10~^, and DOP853 with S = 10"^°, and (b) TM method with r = 0.1, TIDES with 5 = 10"^°, and DOP853 with 
5 = 10-". 



the time evolution of GALIs for a regular 7g orbit with initial condition qi = 0.1, pi = 0, 1 < i < 12, 
for the system ([1]) with = 12, over t = 10^ time units. Further information on these runs are given in 
Table El 



Table 2. Table similar to Table [T] containing information on the performance of the dif- 
ferent numerical methods used for the computation of all the GALIs of the Tq^ orbit with 
initial condition qi = 0.1, Pi = 0, 1 < i < 12, of system ([ij with N = 12, over t = 10^ time 
units (see also Fig. [3]). 
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From Table [2] and Fig. [3] it is seen that for the non-symplectic integrators a one-step accuracy of 
6 = IQ-^'^ is not sufficient for a correct identification of the regular orbit. Results obtained both by 
the TIDES and the DOP853 integrators show a deviation from the theoretically predicted behaviors ([5]) 
after t = 10^ (see the first and third columns of Fig. [3]). To obtain the correct behavior of GALIs over the 
whole integration interval it is necessary to decrease the one-step accuracy to 6 = 10"^^ (see the second 
and fourth columns of Fig. [3]). We note that this decrease in 5 naturally results in a significant increase 
in CPU time. In contrast, the GALIs obtained via the TM method with r = 0.5 require significantly less 
CPU time, and show the theoretically expected behaviors up to t = 10^, although the relative energy error 
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Fig. 3. Time evolution of GALIs for a regular Tq"^ orbit with initial condition = 0.1, Pi = 0, 1 < i < 12 of system ^ with 
A'^ = 12, as computed using non-symplectic schemes. The results are given as colored curves, while the TM method results 
with SABA2C and r = 0.1 are given in grey-scale in the background as reference. 



is \AH/H\ !^ 10"^. 

While the error in energy conservation \AH/H\ grows with time for non-symplectic integrators, it 
remains bounded for symplectic integrators, as can exemplarily be seen in Fig. HI Thus, in case it is 
necessary to integrate beyond t = 10* , one has to further decrease 6 for the non-symplectic methods, in 
order to achieve the same final \AH/H\. This, of course, would result in a further increase in the ratio 
values of the CPU time required by the non-symplectic methods as compared to the TM technique. 



4.2. Searching for motion on low- dimensional tori 

One of the advantages of the GALI method is its capability to identify motion on low-dimensional tori. For 
a regular orbit the largest order k of its GALIs that eventually remains constant deter mines the dimension 



of the torus on whi ch the motion occurs (see Eq. ([5])). This ability was verified in [Skokos et all 12008 



Bountis et all l2009l | . where some particular orbits on low-dimensional tori were considered. 

Since the TM method provides reliable evaluations of the GALIs even for relatively large integration 
steps, and therefore requires little CPU time, we applied this technique to perform a more global investi- 
gation of the FPU-/3 lattice, aiming to trace the location of low-dimensional tori. In particular, we consider 
the Hamiltonian system ([T]) with = 4, for which regular motion can occur on an s-dimensional torus 
with s = 2, 3, 4. According to Eq. ([5]) the corresponding GALIs of order k < s will be constant, while the 
remaining ones will tend to zero following particular power laws. Thus, in order to locate low-dimensional 
tori we compute the GALI^, /c = 2,3,4 in the subspace ((73,^4) of the system's phase space, considering 
orbits with initial conditions qi = q2 = 0.1, pi = p2 = Ps = 0, while pi is computed to keep the total 
energy H constant at H = 0.010075. 

Since the constant final values of GALI^, k = 2, . . . , s, decrease with increasing order k (see for example 
the GALIs with 2 < A; < 6 in Fig. [3]), we chose to 'normalize' the values of GALI^, /c = 2,3,4 of each 
individual initial condition, by dividing them by the largest GALI^ value, max(GALIfe), obtained from all 
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Fig. 4. Time evolution of the absolute value of the relative error in energy conservation for the 7g orbit with initial condition 
Qi = 0.1, Pi = 0, 1 < z < 12, of Hamiltonian ((T| with = 12. The symplectic routine SABA2C uses a step size of r = 0.1, 
while the non-symplectic methods required a one-step accuracy oi S = 10~^^ (see also Table[2|. 



studied orbits. Thus, in Fig. Owe colored each initial condition according to its 'normalized GALI^ value 
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Fig. 5. Regions of different g^. (|13p values, fc = 2, 3, 4, on the (53, 54) plane of the Hamiltonian system ([1} with N = A. 
Each initial condition is integrated by the TM method with SABA2C and r — 0.5 up to t = 10^, and colored according to its 
final (a) g2, (b) 173, and (c) 34 value, while white regions correspond to forbidden initial conditions. Three particular initial 
conditions of regular T2 , T3 and T4 orbits are marked by a triangle, a square and a circle respectively. 

In each panel of Fig. O large values (colored in yellow or in light red) correspond to initial condi- 
tions whose GALIfc eventually stabilizes to constant, non-zero values. On the other hand, darker regions 
correspond to small values, which result from power law decays of GALIs. 

Consequently, motion on 2-dimensional tori, which corresponds to large final GALI2 values and small 
final GALI3 and GALI4 values, should be located in areas of the phase space colored in yellow or light red 
in Fig. [5l[a), and in black in Figs. [SJb) and (c). A region of the phase space with these characteristics is 
for example located in the upper border of the colored areas of Fig. \5\ A particular initial condition with 
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^3 = 0.106, = 0.0996 in this region is denoted by a triangle in Fig. [H This orbit is indeed a T2 regular 
orbit as we see from the evolution of its GALIs shown in Fig. [6|^a). In a similar way, a orbit should be 
located in regions colored in black only in the 34 plot of Fig. Oj^c). An orbit of this type is the one with 
53 = 0.085109, 54 = 0.054 denoted by a square in Fig. [5l which actually evolves on a 3-dimensional torus, as 
only its GALI2 and GALI3 remain constant (Fig. [6jb)). Finally, the orbit with 53 = 0.025, (74 = (denoted 
by a circle in Fig. [5]) inside a region of the phase space colored in yellow or light red in all panels of Fig. [5l is 
a regular orbit on a 4-dimensional torus and its GALI2, GALI3 and GALI4 remain constant (Fig. [6]^c)). It 
is worth noting, that chaotic motion would lead to very small GALI^, and values, since all GALIs would 
tend to zero exponentially, and consequently would correspond to regions colored in black in all panels of 
Fig. [5l Thus, chaotic motion can be easily distinguished from regular motion on low-dimensional tori. 




Fig. 6. The time evolution of GALIs for regular orbits lying on a (a) 2-dimensional torus, (b) 3-dimensionaI torus, and (c) 
4-dimensional torus of the Hamiltonian system ([1} with = 4. The initial conditions of these orbits are marked respectively 
by a triangle, a square and a circle in Fig. [S] 

From the results of Figs. [5] and [6] it becomes evident that the comparison of color plots of 'normalized 
GALIfc' values can facilitate the tracing of low-dimensional tori. The construction of such plots becomes 
a very demanding computational task, especially for high-dimensional systems. Thus, the application of 
the TM method for obtaining such results becomes imperative, since the required computations can be 
performed very efficiently by this method. 

5. Summary and Conclusions 

We compared different numerical techniques for the integration of variational equations of multi- 
dimensional Hamiltonian systems. In particular, we considered the TM method, which uses symplectic 
integrators for the realization of this task, as well as non-symplectic algorithms, like the general-purpose 
Runge-Kutta integrator DOP853, and the TIDES algorithm, which relies on Taylor series expansion tech- 
niques. These methods were applied to the A^D Hamiltonian system ([1]), the FPU-/? lattice, with A^ varying 
from A^ = 4 to A^ = 20. 

We used the numerically obtained solutions of the variational equations for the computation of the 
GALI chaos indicators. The accurate reproduction of theoretically known behaviors of GALIs was used as 
a measure of reliability of the numerical techniques tested. In addition, the CPU time required by each 
method in order to achieve accurate results, was taken into account for the characterization of the efficiency 
of these algorithms. 

The TM method exhibited the best numerical performance in all our simulations, both in accuracy 
and speed. More specifically, we found that the ratio of the CPU time required by the TIDES and DOP853 
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algorithms, with respect to the TM method, for correctly characterizing the nature of orbits, increased with 
increasing N (Fig. [2]). Thus, the TM method should be preferred over the other two techniques, especially 
for studies of multi-dimensional systems. 

A feature of the TM method, which is of significant practical importance, is that it succeeds in finding 
the correct GALI behavior, and consequently determines the nature of orbits correctly, even when large 
integration steps are used, despite the fact that in these cases the energy accuracy is rather low. Therefore, 
the application of the TM method allows the efficient investigation of the dynamical properties of a large 
number of initial conditions in feasible CPU times. As an example, we showed in Sect. 14.21 how the TM 
method can exploit the properties of GALI to efficiently find the location of low dimensional tori in the 
phase space. Possible applications of this approach could be the tracing of quasiperiodic motion in multi- 
dimensional systems, when only a few degrees of freedom are excited. 
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